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Abstract. We discuss the distribution of various estimators for extracting the diffusion constant 
of single Brownian trajectories obtained by fitting the squared displacement of the trajectory. The 
analysis of the problem can be framed in terms of quadratic functionals of Brownian motion that 
correspond to the Euclidean path integral for simple Harmonic oscillators with time dependent 
frequencies. Explicit analytical results are given for the distribution of the diffusion constant 
estimator in a number of cases and our results are confirmed by numerical simulations. 



1. Introduction 



The tracking of single particles is a powerful tool to probe physical and biological processes at the 
level of one macromolecule. In particular, the accumulation of experimental data in recent years 
has allowed to test models of diffusive transport in cells jll ||. Within aqueous compartments, e.g. 
the cell cytoplasm, Brownian diffusion is the basic transport mechanism for proteins pj. Other 
studies, however, have reported subdiffusive behavior both in membranes |l| and in the cytoplasm 
[§), although the microscopic origin of anomalous diffusion remains unclear in this context. Crowded 
environments of the cell may cause slower diffusion than in pure water or other solvents, although 
not necessarily subdiffusion ||. 

Conflicting results have generated a debate on the methodology for determining diffusion laws from 
single particle data, even for simple diffusion |^]. In experiments, trajectories of high temporal and 
spatial resolution are often obtained at the expanse of statistical sample size. Trajectories may be 
few and short due to observation windows limited in space, a rapid decay of fluorescent markers or 
particle denaturation Q . These limitations complicate the determination of the nature of diffusion, 
i.e. a precise estimate of the diffusion constant or an anomalous exponent. 
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In any case, time averaged quantities associated to a trajectory may be subjected to large 
fluctuations among trajectories. In the continuous-time random walk model of subdiffusivc motion, 
time-averages of particle's observables generally are random variables distinct from their ensemble 
averages 0. For instance, the square displacement (after a time lag t) time-averaged along a 
given trajectory differs from the ensemble average JsJ . By analyzing time-average displacements 
of a particular realization, sub-diffusive motion can actually look normal, although with strongly 
differing diffusion constants from one trajectory to an other Q. The Brownian case is different, but 
not as straightforward as often thought. Ergodicity, namely, the equivalence of time and ensemble- 
averages of the square displacement, only holds in this case in the infinite sample size limit. In 
practice, standard fitting procedures applied to finite (although long) trajectories of a same particle 
unavoidably lead to fluctuating estimates of the diffusion constant. Indeed, variations by orders 
of magnitude have been observed in experiments and simple random walk simulations ||. To 
our knowledge, no analytical results are available on the properties of these diffusion constant 
distributions. 

In this article, we present analytical and numerical results on the distributions of the diffusion 
constants estimated from single trajectories. We consider a standard fitting method based on 
time-averaged square displacements as well as other similar procedures amenable to analytical 
calculations. Generally we show that the problem consists of finding the distribution of a quadratic 
functional of Brownian motion with a time dependent measure. 

The first studies of the quadratic functionals of Brownian motion date back to a classic paper of 
Cameron and Martin in 1945 |10| a nd the problem has received much interest in the probability 
community ever since [ fLT[ |l2| , |13|, |14|, |l5| . The formulation of path integrals for quantum mechanics 
provided a powerful tool to analyze this set of problems using methods more familiar to physicists 
[ fl6| |l7j, here the problem appears as the computation of the partition function of a quantum- 
harmonic oscillator with time dependent frequency. Various quadratic functionals of Brownian 
motion have been intensely studied by physicists using a variety of methods. They arise in a 
plethora of physical contexts, for polymers in elongational flows |L9|, a variety of problems related to 
Casimir/van der Waals interactions and general fluctuation induced interactions [^0[ |2l], ^2[ |2^, pi[ , 
where, in harmonic oscillator language, both the frequency and mass depend on time. Quadratic 
functionals of Brownian motion also arise in the theory of electrolytes when one computes the one- 
loop or fluctuation corrections to the mean field Poisson-Boltzmann theory j2f| |26|, ^7], . Finally 
we mention that functionals of Brownian motion also turn out to have applications in computer 
science p9[ . 

In this paper we use the Feynman-Kac theorem to show that the generating function, or Laplace 
transform, of the probability density function of the estimators for diffusion coefficients can be 
expressed as a solution to an imaginary time Schrodinger equation. This Schrodinger equation 
describes a particle in a quadratic potential, whose frequency is time dependent. For the choices 
of time dependent frequency arising in the problem of estimated diffusion constants the resulting 
Schrodinger equation can be solved exactly. The inversion of the resulting Laplace transform to 
obtain the full distribution cannot be carried out exactly, however we are able to analyze the 
behavior of the distribution in both the lower and upper tails, thus giving a rather complete 
analytical description of its behavior. 

In general we find that the main characteristics of the distribution of the estimated diffusion 
coefficient depend little on the fitting procedure used and in all cases its most probable value 
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is much smaller than the correct (average) diffusion constant. The probability of measuring a 
diffusion constant lower than average is actually larger than 1/2 (close to 2/3). 



2. Fits for the diffusion constant of a single trajectory 



Consider a one-dimensional Brownian process B t of variance (B 2 ) — 2D$t = dot. Without 
restricting generality, we set ao = 1 and < t < 1 in the following. If a particular trajectory 
B t is available but ao not known a priori, an estimate a of this parameter can be obtained by 
performing a fit to the diffusion law. Several fitting procedure have been discussed in the context 
of molecule tracking within cells Below, we consider 4 of them. 

One of the simplest method consists in calculating a least squares estimate based on the 
minimization of the sum 



F= f\B?-l(t)] 2 dt, 
Jo 



(1) 



where the diffusion law l(t) can be taken as linear, 

l{t) = a L t, (2) 

or affine, 

l(t) = a A t + b A , (3) 

typically. Given J5 t , the minimization of (0) with respect to the constant(s) yields the least squares 
estimate 

a L = 3 [ tB\dt (FITl), (4) 

for the linear fit, and 

a A = 6 I (2t-l)B?dt (FIT2) (5) 

b A = - 2 I {3t- 2)B 2 t dt (6) 
Jo 

for the affine one. 

Another related method, often used in particle tracking experiments Q and numerical studies [0, 
consists in least-squares fitting the time-averaged square displacement, 5 2 t- For a finite trajectory, 
this quantity is defined as 

8*t= YZ~ t I \Bt + s-B s fd s . (7) 







Due to the ergodicity of normal diffusion processes, at times short compared to 1 the above average 
coincides with the ensemble average (B 2 ) M, i.e., S 2 t — t as t — > 0. However, due to practical 
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Figure 1. Left panel: Square position of a Brownian motion with ao = 1 as a function of 
time and the corresponding diffusion laws obtained with the fitting methods 1, 2 and 4. For this 
example, = 0.318, cia = 0.397 and cimle = 0.338, three values significantly smaller than 
unity. Right panel: time-average displacement calculated for the same trajectory, where Fit 3 
gives a^ 5 ' = 0.274. Only at very short times S 2 t follows the ensemble average dot. The trajectory 
is a random walk of N = 50, 000 steps, with positions x n = X/2=i wnere 1 < n < N and 
li = ±1. In the notation of the text, n/N -> t and x^/N -» Bf. 



limitations, experimental trajectories often have a small number of positions and S 2 t is analyzed for 
all (or a large fraction) of the available intervals t, like in rcf. ||. Similarly, we do not restrict here 
to t <C 1 but fit over the whole time domain < t < 1 instead. As shown by the numerical example 
of Figure (Q-right) for a random walk with N = 50, 000 positions, the expected small t behavior of 
S 2 t can be restricted to a very small interval compared to the total walk duration. Substituting B 2 
by S 2 t in Eq.(Q) and adopting the linear fit, the new estimate simply reads: 

ri 

af = 3 / t 6 2 t dt (FIT3). (8) 
Jo 

Yet another fitting method consists in maximizing the unconditional probability of observing 
the whole trajectory B t , assuming that it is drawn from a Brownian process with mean-square 
displacement at. Namely, the maximum likelihood estimate (MLE), denoted as aMLE, is the value 
of a that maximizes the likelihood of B tl defined as: 



L = f[P a (B u t) = H^ai)- 1 / 2 exp (-^-) 

t=o t=o \ a ' 



(9) 



where P a (x,t) is the probability density of the Brownian process with constant a. By equating 
9 In L/ da to zero, one obtains 

f 1 B 2 

a MLE = dt-i- (FIT4). (10) 
Jo 1 

The estimates given by the four methods above are represented in an example, see Figure (|l|). The 
numerical values are comparable but can differ significantly from unity. 
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Figure 2. Left panel: Distributions of the parameters X = ai, (« symbol) and ctMLE (P 
symbol). tnset: zoom of the same plot, where the solid lines represent the analytical expressions 
(ktj) and (Wq) valid for small X. Right panel: Distributions of the parameters (• symbol) 
along with (solid line) and a a (o symbol). Except for these results are obtained by 

averaging over 2 10 5 random walks with N = 5 10 5 steps. 



3. Numerical results 

The numerical distributions of the random variables at,, clmle, o,^ an< ^ a A are displayed in Figure 

The distributions are highly asymmetric and peaked near X = 0, far from the average value (X) = 1. 
The most probable X is a small positive number in each case, see Table 1. Although estimates of 
X ~ 10 can be sometimes observed, the median of the distribution is lower than (X). Namely, the 
probability of measuring a diffusion constant lower that the correct value is not 1/2, but close to 
2/3 in all four cases. The probability of measuring a negative a a is not zero in the affine method 
(as already noticed in ref. ||) but close to 0.175. Table 1 summarizes the main properties of the 
distribution functions. 

Importantly, and practically obey the same distribution (Figure (Q-right)) , which is 
somewhat unexpected as 6 2 t is a much smoother function than B\. Thanks to this similitude, 
the analytical study of the simpler functional (^|), exposed in the next Section, brings many insights 
on the behavior of . Distributions similar to ours for a£ were determined in ref. || , both 
numerically from random walk simulations and experimentally using R-phycoerythrin proteins in 
mammalian cells. 
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X 




(6) 
a L 


O.MLE 


a A 


{X) 


1 


1 


1 


1 


most probable X 


0.11 


0.16 


0.25 - 0.3 


0.01 


median 


0.54 


0.56 


0.66 


0.42 


lower 5% 


0.086 


0.12 


0.17 


-0.20 


upper 5% 


3.43 


3.33 


2.97 


4.08 


Prob[X < (X)} 


0.683 


0.681 


0.668 


0.683 



Table 1: Main properties of the diffusion constant distributions. 



4. Feynman-Kac formalism for the generating function 



In general the estimated fit parameters discussed above (FIT1, 2 and 4) are quadratic functionals 
of Brownian motion and take the form 



X 



f 

Jo 



B s w(s)ds. 



(11) 



When w(s) > on [0,1] the quadratic functional is positive and its generating function of X, is 
defined by 



/>oo 

= / p(x) cxp(— ax)dx = E [cxp(— crX)] 
Jo 



(12) 



where p(x) is the probability density function of X. In order to compute G we consider the following 
average of a quadratic functional of Brownian motion: 

t-i 



cxp( — g B^w(s)ds) 



(13) 



where the expectation above is for a Brownian motion starting at x at time t. Clearly in this 
notation we have G(a) = \P(0, 0). We now write a Feynman-Kac type formula for ty(x,t) by 
considering how the functional evolves in the the time interval (t,t + dt). During this interval the 
Brownian motion moves from x to x + dB t , where dB t is an infinitesimal Brownian increment such 
that (dB t ) = and (dB^) — dt. Taking into account this evolution we can write to order dt 



*(&,*) = (E- 



exp(— a I Bgw(s)ds) 
Jt+dt 



(1 - dtaw(t)x 2 )) 



(14) 



where the brackets on the right hand side denote the average over dB t . The above may now be 
written as 



V{x,t) = (V(x + dB u t + dt)(l - dtaw(t)x 2 )). 



(15) 



On the distribution of estimators of diffusion constants for Brownian motion 7 

Expanding to second order in dB t and dt, taking the average over dB t and equating the terms of 
0(1) and 0(dt) we obtain 

flT = -2ft? + 0W ^*' (16) 

which looks like a Schrodinger equation in a harmonic, time-dependent potential. The boundary 
condition for this equation is given by Vt'(x, 1) = 1 for all x. 

It is easy to see that the solution of equation ( |l6| ) is given by 

*(x,t) = f( t )exp(-±g(t)x 2 ) (17) 

where 

d i =g2 - 2 ™' . . . . (19) 

with the boundary conditions g(l) = and /(l) = 1. Now we can eliminate the nonlinearity in the 
second equation by setting g = —dh/dt/h which gives 

, df 1 „dh , 

h i + 2 f dt=° (20) 

d 2 h 

- 2awh = 0, (21) 

with the boundary conditions h(l) = 1 and dh/dt(t = 1) = 0. In terms of these functions the 
Laplace transform is now given by G(a) = /(0) = 1/ y/h(Q). We now make a change of time 
variable writing 

dr 



— = y/2w(t)a, (22) 

assuming for the moment that w(t) is positive. In terms of this new temporal variable equation 
( pi] ) can now be written as 

^ + JlL^_ /l = o (23) 
dr 2 (jl) 2 dr 

In the class of problems we study in this paper (see Eqs.([|), and (|lo|)) the form of w is 

w(t) = (At + C) a , (24) 
with A and C two constants. From this we can choose r to be 



\A\(a + 2) 

and equation (^3|) becomes 

d 2 h a dh 

dr 2 + (a + 2)t dr 



Vr (At + C)^ (25) 



h = 0. (26) 
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The general solution to this equation can be shown to be 

h(r) = (DK^_ (r) + EI^_ (r)) , (27) 

where K v and I v arc modified Bessel functions J3(i| . The coefficients D and E are determined from 
the boundary conditions h(ri) — 1 and dh/dr = at t\ = r(l) = \/8a(A + C) 2 ^ /\A\(a + 2). 
Solving for D and E and using standard identities for Bessel functions |$0| we find that at 
t = r(0) = y/j&C 2 * 1 /\A\(a + 2) 

Hr ) = t~*t?& (/^(n)^ Wr ) + (n)J , (r )) , (28) 

and thus 

(29) 



G(a) = 



5+1 

i + 2 



n)X 



(to) + ^_h±i (ti) (r ) 

. + 2 ^+2 V 7 a + 2 



5. Asymptotic analysis for the probability density function 



The general result equation (|29|) simplifies in the case where To = 0, i.e. when C = 0, which is 
the case for FIT1 (linear) and FIT4 (MLE). In this case the probability density function of the 
estimator of the diffusion coefficient p(x) has support on [0, oo). We start by analyzing the behavior 
of p(x) at small x. 



We proceed by using the small argument expansion of K v for v > 0: 

K v {z)^\t{u){\z)- v 
to obtain the exact result 



G(a) 



r(— ^) 



\/2crA a 
a + 2 



i+i 

a + 2 



/ \/8o\A Q 
Si \ a + 2 



(30) 



(31) 



The moments of X can then be extracted using the scries expansion for modified Bessel functions 
[p0| which gives 



G(a) = 



( 2oA*\ ' 



(32) 



Without loss of generality we set A = 1 and find the first two moments of X to be given by 

W = zr^ ( 33 ) 



and thus 



(X 
(X 2 



a + 2 
, 3a + 7 



(a + 2) 2 (a + 3) 
2 

= (a + 2)(a + 3) 



(34) 
(35) 
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Figure 3. Laplace transforms of the distributions of X = and a a,r r n (cases {A = 3, a = 1} 
and {A = 1, a = —1}, respectively). The solid lines are given by Eq.felh; the points represent 
the simulations results. 



In FIT1 and FIT4, a single estimator for the diffusion constant has the form 

X a = {a + 2)X = {a + 2) [ t a B?dt, (36) 

Jo 

with a = 1 and —1, respectively, which gives 

<^) C = 2(1--L-), (37) 
a + 6 

From this we see that the MLE estimate of the diffusion coefficient has a variance (^ 2 _ 1 ) = 1 where 
as the simple linear fit has a larger variance (Xf ) = 3/2. Of course these variances can be computed 
directly and the above analysis serves as a check on our formalism to compute the full probability 
density function. 

An interesting comparison can be made with the estimator X ep which uses just the final value of 
the mean squared displacement 

X ep = Bl, (38) 

here we find the variance 

(X 2 ep )c = 2, (39) 

which is clearly bigger than all the integral estimators above. Before embarking on inversion of 
the generating function G(tr) to obtain the probability density function p(x), a simple check of 
our results is to numerically compute G(o~) from our simulation data. In Figure (0) are shown the 



Laplace transforms G(a) obtained from both Eq.(31 ) [or (|32j)] and the numerical distributions p(x) 
we see that the agreement is perfect. 



The behavior of X at small values (when it is always positive) can be extracted by examining the 
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characteristic function, or equivalently the Laplace transform of the probability density function 
p{x) of X. Using the large z asymptotic expansion 

^Lexp(z) (40) 



/2ttz 

and setting A — 1, we find for large a 



The behavior of p(x) at small x can now be extracted by noticing that the integral 

f°° d 

I = exp(— ax) exp( )x c dx 

Jo x 

is dominated by its value at small x and thus can be evaluated by the saddle point method as 



(42) 



2c+l 

/- ^exp(-2Vad) (^j 4 (43) 
from which we deduce that for small x 

p(x) ~ 7r-ir-i(-l^)(a + 2)-M ^-Ig^ exp ^_ __L_ ) . (44) 

From this we obtain the probability density of X = X a [Eq.(|36|)] at small x to be: 

p a (x) ~7r-3r-5(^_)( a + 2)-3fS% ^-jg^expf- ) . (45) 

a + 2 \ 2(a + 2)x / 

The distribution exhibits an essential singularity at x — 0, as expected from the general asymptotic 
result of Shi |n|. For the linear fit estimate (a — 1), Eq.(|45|) gives 

Pi{x) ~ ci x~ 12 exp ( — — ) (46) 



G.i: 



Cl = 3-A7r-3r(i)-5 ~ 0.29035..., (47) 



with 

v 3 



and for the MLE (a = -1) 

~ c_i aT^ exp f _ ^J ( 48 ) 

with 

c_i = 7T"3 w 0.75112... (49) 

The expressions above compare well with the simulation results at small x (Figure (^-left), 
inset). The distributions ( |46] ) and ( |48| ) actually present a maximum at a;* = 2/17 « 0.118 and 
a;* = 2/7 ps 0.286, respectively. Despite that the asymptotic results start to fail when x becomes 
too large, these values are still in good agreement with the most probable values of Table 1. A 
more detailed comparison in the small x regime is displayed in Figure (^), where p(x)x@ obtained 
from the numerics is plotted as a function of 1/x, with j3 = 17/12 and 7/4. The behaviors at 
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0.1 



0.01 



0.001 



1e-04 



1e-05 




Figure 4. Rescaled numerical distributions p(x)x^ with = 17/12 (linear fit, black dots) and 
= 7/4 (MLE fit, diamonds) as a function of 1 At. The solid lines are the analytical forms 



ci exp(— J-) and c_iexp(— -!-) from Eqs.(k 



and 



), respectively. 



large arguments are nearly indistinguishable from the exponential laws predicted by Eqs.([46|) and 
©■ 

In order to extract the behavior of the probability distribution for large x we need to examine the 
singularities of the generating function G(a) for a < 0, in this regime 

-1/2 



G(a) = 



■i±i 

> + 2 



a + 2' \ a + 2 



J 



\ a + 2 



(50) 



from the identity J„{z) = YX=a(- 1 ) k ( z l 2 Y k+V /W-T{k + v + 1)]. This Bessel function of the first 
kind oscillates and has simple zeros, at these zeros G diverges. Let us denote u* as the lowest 
positive zero of J «+i (u). When u = y/ '8\a\A a / (a + 2) -+ u* from below, 



J a+l (u) 



-1/2 



\ 



21(7* 



u*\J'_ a+1 (u*)\ 



where a — > a* = —u* 2 (oe + 2) 2 /(8A a ) from above. We now note that 

exp(— lux) 



dx 



■ exp(o~x) 



sjx y oj — G 

for to > a. Comparing Eqs. (J5l|) and (B2J) , one deduces from (KQ) the large x behavior: 



p(x) 



\a-\2 e -\cr*\x 



(51) 



(52) 



(53) 
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Figure 5. Left panel: Rcscaled numerical distributions p(x)x 1 i 12 for the linear (black dots) and 
MLE (diamonds) fits. The solid lines are the exponential laws from Eqs.(^) and (^), Right 
panel: Same quantity for the affinc fit (FIT2). The solid lines are the asymptotic forms at large 
x and large —x, see Eqs.(Ej) and (Ej). 



For the linear fit (a = 1, A = 3), one finds u* = 1.2430... and 

-0.5794k 

p x {x) « 1.1675 — (54) 
V 2nx 

whereas for the MLE (a = -1, A = 1), u* = 2.4048... and 

-0.7228a; 

p_i (x) « 1.5212 — (55) 
V 2nx 

These asymptotic expressions are compared with the numerical results in Figure (|5]deft). 

The interpretation of this result is rather straight forward, if we consider a Gaussian random Y 
variable of mean zero and variance j 2 then the probability distribution is 

^^T^? 6 ^"^- (56) 

Now defining Z — Y 2 we find that the the probability density function of Z is 

exp(-^r) 

pz{x) = -^r- (57) 

which has the same functional form as equation (|53|) . This means that for large values of x the 
random variable X has the same distribution as a squared Gaussian random variable. This is 
not surprising as the variable X can be viewed as an infinite sum of Gaussian random variables. 
Note that the full probability density function for the end point estimator X ep is given by (as 
7 2 = 2) 

Pep (x) = eXp( ^ } , (58) 
and so the distribution of this simple estimator decays much more slowly that the two integral 
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estimators discussed above. 

In the case of the affine fit, FIT2, both the estimators aA and bA, defined in equations (||) and (||), 
can be negative as the respective functions w change sign. The probability density function is thus 
two sided. When tq becomes imaginary in Eq. (|28|) , this solution must be modified by substituting 
h/3( T o) and /_i/ 3 (t ) by -Ji/ 3 (|t |) and J-i/ 3 (\t \), respectively @. In turn, when t\ becomes 
imaginary, J 2 /3(n) and I_ 2 / 3 (n) are replaced by J 2 / 3 (|n|) and J_ 2 / 3 (|ti|), respectively. For large 
x > the probability density function can be obtained from the closest zero of h(o~) from zero in 
the negative direction, denoted by a*_, and the analysis above goes through to give 

e -l CT - \ x 

p(x) » , with |ct* | = 0.4596... and A- = 0.9239..., (59) 

\/2ttx 

in the case X = a A . For the variable X = b A , one finds |cr* | = 3.2229... and A_ = 1.4734... As 
X can become negative we also have zeros of h(a) for positive values of a; now if the first of these 
zeros from the origin is then the same analysis as above implies, for x < 0: 

p(x) » A , - + with ct* = 4.2439... and A + = 0.8381..., (60) 

v 2?r M 

in the case X = cia- For X = b A , one obtains a*^ = 2.4485... and A + = 1.1886... These asymptotic 
results are tested in Figure (fright) on the two sided distribution arising for both the coefficients 
a a and b A , showing very good agreement. 

6. Conclusion 

We have shown that a general class of statistical estimators that can be used to extract diffusion 
constants from the squared displacement of single Brownian trajectories are in fact quadratic 
functionals of Brownian motion. Numerically we have seen that such estimators have a tendency 
to yield values which are typically lower than the correct average value. In addition we have seen 
that the statistics of the estimated diffusion constants from these trajectories resemble closely those 
obtained from fitting the time averaged squared displacement 6 2 t, defined in equation (0), despite 
the fact that the resulting trajectory appears much more regular than an unaveraged Brownian 
squared displacement, as demonstrated in Figure (JjJ-right). An interesting and outstanding problem 
would be to carry out our analysis for estimators of type Sf. Such an extension is clearly desirable 
as it deals with a quantity more commonly used in single particle tracking experiments. However 
from a technical point of view the resulting path integrals, while being for quadratic functionals of 
Gaussian processes, are highly non-local in time and it is probable that their evaluation will require 
the introduction of new mathematical methods. 

Our final analysis was only limited by the problem of carrying out a full Laplace inversion of the 
generating function G(a) to obtain the full probability density function. However we point out 
that the generating function is actually easy to estimate from numerical data for the purpose of 
comparison with our analytical results, as demonstrated in Figure ([3]). In addition the generating 
function can be inverted analytically in certain asymptotic regimes. When the estimator is always 
positive, and consequently p{x) = for x < 0, the behavior of p(x) for small x can be extracted. 
We find that it has an essential singularity at x = and a maximum value, this estimate of the 
maximum value is in good agreement for the most likely value of x coming from the full probability 
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density function. For positive estimators the large x behavior of p(x) turns out to be that of a 
squared Gaussian random variable, reflecting that fact that the estimator itself is an infinite sum 
of Gaussian random variables. This remains true when the estimator can have negative values, i.e 
when w(t) can change sign. In this case the probability density function for X is that of a Gaussian 
squared for large x as is that of —X for large negative x. 

Finally new methods are being introduced into single particle trajectory analysis to estimate 
diffusion constants and exponents associated with anomalous diffusion, for instance methods based 
on the mean maximal excursion and it would be interesting to examine the distributions 
associated with such estimators. 
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